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SUMMARY 


Applying a finite difference approximation to a biharmonic equation results in a very 
ill-conditioned system of equations. This paper examines the conjugate gradient method used in 
conjunction with the generalized and approximate polynomial preconditionings for solving such linear 
systems. An approximate polynomial preconditioning is introduced, and is shown to be more efficient 
than the generalized polynomial preconditionings. This new technique provides a simple but effective 
preconditioning polynomial, which is based on another coefficient matrix rather than the original 
matrix operator as commonly used. 

INTRODUCTION 

The conjugate gradient (CG) method used in conjunction with a suitable preconditioning has 
been widely accepted as one of the most efficient iterative procedures for solving large and sparse 
systems of linear equations. The basic CG algorithm can be efficiently implemented on a vector 
computer. However, the traditional preconditioning technique, which is generally based on an 
incomplete LU factorization procedure (Meijerink and van der Vorst, 1977) is not suitable for coding 
on a vector computer. The solution for the preconditioning step is obtained via forward and backward 
substitutions, which are difficult to vectorize because of the recursive nature. To overcome this 
problem, the use of a matrix polynomial as a preconditioner has been proposed (Dubois et al., 1979). 
The preconditioning operator is taken to be an approximation of the matrix inverse, which can then 
be obtained by taking sufficient number of terms in the matrix polynomial expansion. The main 
advantage of such an approach is that the preconditioning step can be implemented essentially via 
matrix by vector multiplications, which do not require forward and backward substitutions. 
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A very important practical consideration in the preconditioned CG algorithm is that the 
preconditioning step should be easily computed. Notice that the computational work for each basic 
CG iteration essentially comprises of one matrix by vector multiplication, two vector inner products, 
and three vector additions. For the incomplete LU preconditioning procedure, the amount of work 
due to the preconditioning is about one matrix by vector multiplication. However, it becomes k matrix 
by vector multiplications when k terms are kept in the polynomial preconditioning technique. Hence, 
the computational work per CG iteration when using a polynomial preconditioning may be dominated 
by the cost required for the preconditioning step when k is greater than one. An earlier paper (Wong, 
1987) presented an approximate polynomial preconditioning technique that may be more efficient 
than the traditional polynomial preconditionings. For matrix problems resulting from high-order 
discretization methods applied to elliptic partial differential equations, this method provides a 
substantial saving in computing time compared to that with the usual polynomial preconditioning. The 
basic idea of this technique is that the preconditioning polynomial is derived from another coefficient 
matrix that contains far fewer non -zero elements than the original matrix operator. 

This paper demonstrates how an approximate polynomial preconditioning can be applied to 
matrix equations resulting from a biharmonic problem. Although our discussion will be restricted to 
symmetric and positive definite matrices, the extension to general non -symmetric matrices with 
positive definite symmetric part is straightforward. 


RESULTS AND DISCUSSION 


and 


Consider the Dirichlet problem for the biharmonic equation in a square D with boundary 9D 
A 2 u(x,y) = f(x,y) in D (1) 

u(x,y) = gj(x,y) 


u n (x,y) = g 2 (x,y) 


on 0D 


where A 2 u = u + 2u + u and u denotes the normal derivative of u(x,y). 

aaaa aa yy yyyy n 


When Eq.(l) is discretized by a finite difference scheme with a uniform mesh h, the result is 
a system of linear equations 

Bu = b (2) 

where B is a symmetric matrix with coefficients (except those close to the boundary) given by the 
following 13 -point finite difference stencil 
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It should be pointed out that although B is a symmetric and positive definite matrix, it does 
not possess the property of diagonal dominance. Furthermore, B is very ill conditioned and has a 
condition number proportional to h' 4 . Hence, the application of an iterative procedure, such as the 
CG algorithm, results in a very slow rate of convergence. 

The basic preconditioned CG algorithm is outlined in Section 1, and the approximate 
polynomial preconditioning technique is described in Section 2. The calculation of parameters to 
improve the preconditioning polynomial is presented in Section 3. The choice of a special initial 
solution to accelerate the CG algorithm is discussed in Section 4. From the computational experiments 
given in Section 5, it appears that the use of the CG algorithm in conjunction with an approximate 
polynomial preconditioning results in a highly efficient iterative procedure for the solution of a 
biharmonic equation on a vector computer. 

1. PRECONDITIONED CONJUGATE GRADIENT ALGORITHM 

Let the preconditioning operator M be a non-singular matrix. Eq.(2) can be rewritten as 

BM^Mu = b (4) 

The CG algorithm applied to the preconditioned system (4) is summarized as follows: 
Initialization: 

Step 1. Choose an initial solution u° 

Step 2. Compute residual vector r° = b - Bu° 

Step 3. Set direction vector p° = M'V 3 
Iteration: For n=0,l,2,..., until convergence, do 

Step 4. a n = (r I1 ,M' 1 r n )/(p n ,Bp n ) 

Step 5. u n+1 = u n + a n p n 

Step 6. r n+1 = r n - « n Bp n 

Step 7. /? n = (r n+1 ,M' 1 r n+1 )/(r n ,M’ 1 r n ) 

Step 8. p n+1 = M' 1 r I1+1 + 0 n p n 
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The inner product is defined as (x,y)=x T y for any vector x and y. When M=I, it becomes 
the basic CG algorithm applied to the original system (2). 

The basic CG algorithm is suitable for coding on a vector computer such as the CDC CYBER 
205. Notice that steps 5, 6, and 8 are the well-known linked triad (i.e., vector + constant • vector) 
operations, and the two inner products required in step 4 and 7 can be computed by using the 
Q8SDOT routine available on the CYBER 205 computer. The main computational work in each CG 
iteration is one matrix by vector multiplication Bp, which can be implemented almost entirely by the 
linked triad operations (Madsen et al., 1976). Consequently, the work for each preconditioned CG 
iteration is given by 

W = W b + W p (5) 

where W fe is the work for the basic CG algorithm and 

W b = 2*IP + 3*LT + 1*MV (6) 

Here IP, LT, and MV denote operations for the inner product, linked triad, and matrix by vector 
multiplication. W p is the work resulting from the preconditioning step, and W p =0 if M=I. 

It is well known that the convergence rate of the CG algorithm can be substantially improved 
by the application of a good preconditioning matrix M. The important considerations, in the selection 
for M are that M' 1 should be a good approximation to B' 1 in some sense, and that the 
preconditioning step M'*r can be conveniently computed. A popular choice for M, which can be 
efficiently implemented on a vector computer, is based on a Neumann series for approximating a 
matrix inverse (Dubois et al., 1979). This procedure is generally referred to as a polynomial 
preconditioning technique. 

Suppose the original matrix B has been symmetrically scaled so that b u = 1.0 for all i, B can 
then be expressed as 

B = I - G (7) 

where G is symmetric with g.. = 0.0 for all i. Since B' 1 
truncated matrix series expansion in G; that is, 

k 

M' 1 = I G 1 

i = 0 


= (I-G)' 1 , M' 1 can be obtained via a 

( 8 ) 
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For convergence, the spectral radius of G must be less than one. 

Unfortunately, when B is resulting from a bihannonic equation, it can not be conveniently 
scaled in the form given in (7) while keeping G symmetric. Due to the boundary conditions, b- is not 
constant: it has a value of 20 for all i in the interior points, and it becomes 21 or 22 for i 
corresponding to the boundary points. 

To alleviate this difficulty it is natural to generalize the polynomial preconditioning technique. 
One approach is by considering the matrix splitting of B. Let B=D-Q, where D is the diagonal 
matrix, then 

k 

M 1 = (Z G j ) D' 1 (9) 

i=0 

where G = D -1 Q. Notice that the preconditioning matrix M is symmetric, and the spectral radius of 
G is less the one if and only if D+Q is positive definite (Adams, 1985). Another approach is by 
expressing B is the following form: 

B = o [I-(I-(B/u))] = u(I-G) (10) 

where u is a parameter. It is not hard to verify that if «>Xmax/2 (where Xmax is the largest 
eigenvalue of B), the spectral radius of G is less than one. A simple choice for « is o= | |B 1 1/2, where 
llBlI is taken to be the maximum row sums of the matrix B. Hence, B 1 =w' 1 (I-G)' 1 , and the 
generalized polynomial preconditioning matrix can be defined as 

k k 

M 1 = Z G' = Z (HB/o)) 1 (11) 

i=0 i=0 

The main advantage of this approach is that the original matrix B does not need to be symmetrically 
scaled. Furthermore, it can be applied to non -symmetric problems and matrices which are not 
diagonally dominant. It is easy to verify that u=l for scaled matrices resulting from a finite 
difference approximation to Laplace operators. For these problems, the generalized polynomial 
preconditioning (11) reduces to the simple form given in (8). 

2. APPROXIMATE POLYNOMIAL PRECONDITIONING 

The generalized polynomial preconditionings described in the previous section can be expressed 
as 

. k 

M' 1 = (Z G 1 ) N (12) 

i=0 
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In order to avoid confusion between these two versions, the following abbreviations will be used: 

(1) DPPl where G=D' 1 Q, N=D' 1 , and B=D-Q 

(2) DPP 2 , where G=I-(B/o), N=I, u= ||b|(/2 

Here, DPP denotes the direct polynomial preconditioning (i.e., the polynomial series is based on the 
original matrix). 

Notice that M' 1 becomes a good approximation to the matrix inverse B' 1 as the number of 
terms k increases in the polynomial expansion (12). The preconditioning step M _1 r required in each 
CG iteration can be calculated via 

M -1 r = (I + G + G 2 + ... + G k ) r (13) 

where r =Nr. The computational work due to the preconditioning is then essentially given by 

W = k*LT + k*MV (14) 

p 

if k terms are kept in the preconditioning polynomial. Since G has the same structure as the original 
matrix B, the arithmetic operations for one matrix by vector multiplication Bp required in W b 
(Eq.(6)) is almost the same as that for Gr needed in W p (Eq.(14)). The computational work for 
each CG iteration given in Eq.(5) may then be dominated by W p if k is greater than one. 
Consequently, reduction in the cost for W p will result in an overall improvement for the 
preconditioned CG iterative process. 

An earlier paper (Wong, 1987) shows that for certain matrix problems it is possible to obtain 
an efficient polynomial preconditioning in which the polynomial expansion is based on another 
coefficient matrix that consists of far fewer non-zero elements than the original matrix operator. To 
distinguish from the traditional polynomial preconditioning, this new approach will be referred to as 
an approximate polynomial preconditioning technique. 

Consider B as the matrix resulting from the difference approximation to a biharmonic 
equation. It can then be shown that 

B = L 2 + E (15) 

where L is the difference approximation to the Laplace operator, and E is a matrix with small rank. 
Hence, when solving the biharmonic equation Bu=b, a polynomial preconditioning can be considered 
in which the matrix polynomial is based on L rather than the original matrix B. The approximate 
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( 16 ) 


polynomial preconditioning (APP) can be expressed as 

. k 
M' 1 = Z G 1 
i=0 

where G=I-(L/w). Since L is the Laplace operator, w=l, and G=I-L. The attractive feature in using 
APP instead of DPPi or DPP 2 is that the computational cost for W p is substantially reduced when k 
is large. Notice that the cost for one matrix by vector multiplication Gr based on APP is less than 
half of that of the DPPi or DPP 2 methods. This is due to the fact that for two dimensional problems, 
the matrix B has 13 non -zero diagonals, whereas L has 5 non -zero diagonals. 

3. PARAMETER SELECTION 

The basic generalized and approximate polynomial preconditionings have been described in the 
previous sections. More efficient polynomials, however, can be achieved by introducing appropriate 
parameters into the preconditioning polynomial. Several authors have discussed ways of finding the 
parameters. For more details see Johnson et al. (1983) and Saad (1985). This section briefly outlines 
how the parameters are determined to improve the approximate polynomial preconditioning. 

By introducing the parameter y i into the APP given in (16), 

k 

M' 1 = Z y. G 1 (17) 

i=0 

where G=I-L. The parameters Yj are chosen so that the eigenvalue distribution of the preconditioned 
system M -1 B is more favourable to the CG method. In general, this will also result in a smaller 
condition number for M -1 B. To satisfy these requirements, the choice of Yj should minimize 
| ll-M^B 1 1. However, it is more convenient to consider | |l-M _1 L 2 1 1 by ignoring the contribution of E 
given in Eq.(15). 

Recall that M' 1 defined in Eq.(17) is a polynomial of degree k in G. If M' 1 =P k (G), it is 
then clear that M _1 L 2 is a polynomial of degree k+2; hence, 

M _1 L 2 = Q k+2 (G) = P k (G)(I-G) 2 = Z yp\l-G) 2 (18) 

i=0 

Since the spectral radius of G is less than one, that is, 

-1 < \ < - < \ < 1 (19) 
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where X. denotes an eigenvalue of G. The corresponding eigenvalues of M^L 2 are then given by 

<WV - j o (1-V 2 (2°) 

A good choice for y. is to make the quadratic norm of 

+ 1 

J tQ k+2 ( x ) ‘ 1 f dX (2D 

-i 

as small as possible. Applying a least square method to minimize the integral (21), it leads to a 
normal equation 


Cy = c (22) 

where C is a dense matrix of order k+1. Alternatively, the method of orthogonal polynomial 
expansion can be applied to construct a three term recurrence relation for the preconditioning 
polynomial. The advantage of this approach is that it does not require solving the ill-conditioned 
normal equation. The following are the resulting first three polynomials: 

k=l, P^G) = I + (7/8)G 

k=2, P 2 (G) = (49/40)1 + ( 98/40 )G + (63/40)G 2 

k=3, P 3 (G) = (25/24)1 + (72/24)G + (117/24 )G 2 + (66/24)G 3 

The same procedure can be applied to determine suitable parameters for the DPP 3 and DPP 2 
preconditionings, which will not be described here. 

4. SPECIAL INITIAL SOLUTION 

An initial solution vector u° is required to start the CG iterative procedure. In general, the 
rate of convergence of the CG method is not sensitive to the initial guess u°. It has been found that 
when the preconditioned CG algorithm is applied to the elliptic difference equations with u°=0.0, or 
u° consisting of uniform random numbers, or u°=M' 1 b where M 1 is a preconditioning polynomial 
which approximates the original matrix inverse, there is little difference in the numbers of iterations 
required for convergence (Wong, 1987). 

For the solution of biharmonic equations Bu=b, however, it is possible to determine an initial 
vector u°, which gives an improved convergence rate for the CG method. Recall that Eq.(15) stated 
that B=L 2 +E. Now, ignoring the matrix E, u° can be obtained by 

L 2 u°=b (23) 
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The solution for u° can be found via solving 
Ly = b 

and . (24) 

Lu° = y 

Since L is the Laplacian difference matrix, the solutions for Eq.(24) can be computed very efficiently 
by the CG method with a polynomial preconditioning. The question to be asked now is: Whether this 
choice for u° provides a good starting initial solution for the CG method? 

Suppose u is the exact solution of the biharmonic equations, and let 

u° = u + e (25) 

The initial residual vector is then given by 

r° = b - Bu° = b - B(u+e) = -Be (26) 

Bu° = Bu + Be 

(L 2 +E)u° = b + Be (27) 

Eu° = Be 

r° = -Eu° (28) 

Therefore, u° provides a small initial residual as long as the norm of E is small. However, it should 
be pointed out that this is not the only factor which results in a better rate of convergence for the CG 
method. The initial solution u° also has a special property that leads to improved convergence rates. 
Notice that u° is the solution of L 2 u°=b, and it is well known that the eigenfunctions of the Laplace 
and biharmonic operators are closely related (Birkhoff and Lynch, 1984). Consequently, not only that 
the resulting initial residual vector is small, but it also contains components of the eigenvectors of the 
biharmonic operator. A faster convergence rate can therefore be achieved when the CG method is 
used with this particular choice for u°. This remark is verified by the computational experiments. 

5. COMPUTATIONAL RESULTS 

In order to determine the effectiveness of the approximate polynomial preconditioning, 
consider the solution of a biharmonic equation (1) over a square region with 0<x,y<l. The 
functions g^x.y) and g 2 (x,y) are chosen so that the exact solution is given by 

u(x,y) = x 2 y 2 (x-1.0) 2 (y-1.0) 2 (29) 


From Eq.(25), 


Hence, 
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Eq.(l) is discretized on a uniform mesh h. and the total number of equations is given by EQ=n 2 , 
where h=l/(n+l). Numerical experiments have been performed on the CDC CYBER 205 computer 
with 2 -pipe and 64-bit arithmetic. For the results reported in this section, the CG iterative process is 
terminated when the maximum | r 4 1 is less than 10' 10 , where r is the residual vector. 

5.1 DIRECT VERSUS APPROXIMATE POLYNOMIAL PRECONDITIONINGS 

The performance of the preconditioned CG algorithm with various preconditionings is 
examined for h= 1/100. The initial solution vector u° is composed of random numbers uniformly 
distributed in [0,1]. In Table 1, results for DPP p DPP2 2 , and APP preconditionings are compared, 
where k is the number of terms in the polynomials. Note that no parameter is introduced into the 
polynomial preconditionings. 


Table 1. Number of iterations for EQ=9801 
(No parameter in preconditionings) 


k 

0 

1 

2 

3 

4 

DPP 1 

5512 

2777 

4260 

2683 

4054 

dpp 2 

5492 

2777 

3162 

1954 

2439 

APP 

5492 

2192 

3206 

1690 

2504 


From the results reported in Table 1, it is clear that a large number of iterations is required for 
convergence when the CG method is used without preconditioning (i.e., k=0). It is also observed that 
DPP 2 provides better performance than the DPPj preconditioning. Hence, our attention for direct 
polynomial preconditionings will be focused on DPP 2 , and it will hereafter simply be referred to as 
DPP. 


As previously mentioned, more efficient preconditioning can be obtained by introducing 
appropriate parameters into the polynomial. The effect of the parameters y^ on DPP and APP 
preconditionings is shown in Table 2. The parameters are calculated by the least square method as 
described in Section 3. The corresponding CPU time in seconds is given in Figure 1. 
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Table 2. Number of iterations for EQ=9801 
(Parameterized preconditionings) 


k 


1 

2 

3 

4 

5 

10 

20 

30 

DPP 

5492 

3024 

2239 

1712 

1416 


685 

379 

261 

APP 

5492 

2335 

1448 


778 

623 

299 

142 

97 


Several interesting phenomena can be seen in these results. First, parameterized 
preconditionings provide continuous improvement as the number of terms k increases unlike those 
based on the basic polynomial preconditionings, as is seen in Table 1. By comparing the results 
corresponding to k=5, 10, 20, and 30, it is noted that the number of iterations required for 
convergence is reduced by half as k increases by a factor of two. Secondly, the APP preconditioning 
should be preferred. Not only is the number of iterations for a fixed value of k is smaller for the APP 
than for the DPP method, but the computing time is also significantly smaller. The standard CG 
method with no preconditioning requires 18.5 seconds for convergence, whereas the best results for 
the CG+DPP and CG+APP methods are 13.1 and 2.2 seconds, respectively. 

The effectiveness of the polynomial preconditionings can be best illustrated by examining the 
eigenvalue spectrum of BM' 1 . Figure 2 gives the original spectrum of the biharmonic operator (i.e., 
M=I) for n=16 (i.e., EQ=225). In Figures 3 and 4, we show the corresponding eigenvalue spectrum 
of the preconditioned operator using DPP and APP preconditionings for k=l, 5, 10, and 15. The 
condition numbers C=Xmax/Xmin are also reported in each case, where Xmax and Xmin are the 
largest and smallest eigenvalues. In Figures 2 to 4, each eigenvalue is plotted by a vertical bar. From 
these figures, it is seen that when the number of terms k increases for both DPP and APP 
preconditionings, not only are the condition numbers reduced but the eigenvalues are also more 
clustered and thus more favourably distributed for the CG method. 

It is evident that of the two preconditionings considered here, the APP method is the most 
effective for the biharmonic problems. This observation is somewhat unexpected, because it might be 
thought that the use of the DPP method (which is based on the original matrix operator) would give 
a better eigenvalue spectrum for the CG method, and hence would require a smaller number of 
iterations. However, it may not be as efficient as the APP because of the large amount of work would 
be needed to perform the preconditioning steps. Numerical results given in Figures 3 and 4, on the 
other hand, reveal that for a given number of terms k in the preconditioning polynomials, the 
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condition number for the preconditioned system based on the APP method is actually smaller than 
that based on the DPP method. Consequently, the APP preconditioning results in a faster rate of 
convergence for the CG iteration, as is reported in Table 2. Moreover, since the computational work 
for each preconditioned CG iteration based on the APP method is much smaller than that for the 
DPP preconditioning, a substantial improvement in the performance of the CG method is achieved. 

5.2 EFFECT ON SPECIAL INITIAL SOLUTION 

Table 3 presents the results of the DPP and APP methods for EQ=9801, in which the initial 
solution is obtained via solving L 2 u°=b as described in Section 4. The corresponding CPU time in 
seconds is shown in Figure 5, where the time includes that for finding the initial guess u°. 


Table 3 . Number of iterations for EQ=9801 
(Parameterized preconditionings and special u°) 


k 

0 

1 

2 

3 

4 

5 

10 

20 

30 

DPP 

1379 

758 

561 

430 

356 

303 

174 

95 

66 

APP 

1379 

596 

370 

262 

207 

169 

88 

47 

33 


The results presented in Table 3 and Figure 5 are compared to those of Table 2 and Figure 1, 
which are based on an initial solution where u° consists of random numbers in [0,1]. The results 
clearly illustrate the advantage of using the special choice of u°, namely their ability to substantially 
reduce the number of iterations required for convergence. Consequently, a significant saving in 
computing time is achieved. The best results for solving EQ=9801 for DPP and APP methods are 3.8 
and 1.2 seconds compared to 13.2 and 2.2 as reported earlier. 

5.3 COMPARISON OF CG AND CG + APP METHODS 

The performance of the standard CG without preconditioning and that with APP 
preconditioning together with the special initial starting solution u° is now investigated. Table 4 
compares the number of iterations required for convergence for various numbers of equations. For 
the APP method, k is fixed to be 25 for all cases tested here. The corresponding CPU time is also 
given in Figure 6, where CG+APP is the total time required to solve Bu=b and the time taken to 
compute u°=(L 2 )' 1 b. 
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Table 4. Number of iterations 


n 

100 

150 

200 

250 

EQ 

9801 

22201 

39601 

62001 

CG 


11674 

19568 

31462 

CG+APP 

32 

71 

126 

170 


These results speak for themselves. The standard CG method is not practical for solving 
ill-conditioned biharmonic equations due to their very slow convergence rates. However, significant 
improvement is clearly achieved when it is applied in conjunction with the APP technique and with 
the special choice for u°. For EQ= 62001, the computing time for CG + APP is 28.5 seconds, whereas 
the CG method requires 686.1 seconds; this represents a reduction by a factor of 24. 

6. CONCLUSIONS 

Many practical problems in elasticity and fluid dynamics require the solution of equations 
involving biharmonic operators. Discretizing the equation by a finite difference method results in a 
very ill-conditioned system of linear equations. The computational experiments reported in this paper 
show that the conjugate gradient method with no preconditioning is not practical because of its slow 
convergence rates. However, the convergence rate can be substantially improved by the use of a 
preconditioning technique. Of the two versions of the generalized polynomial preconditionings DPP 1 
and DPP 2 presented in this paper, the DPP 2 method is preferred. This method provides a better 
performance than the DPP X method, and it does not require the matrix to be symmetrically scaled or 
to possess the property of diagonal dominance. 

Another technique introduced, namely, the approximate polynomial preconditioning (APP), 
provides a significant improvement over the traditional polynomial preconditioning. The conclusion 
favouring the APP method is based largely on the grounds of computational work involved in the 
preconditioning step. The application of the CG algorithm in conjunction with the APP technique 
results in a highly efficient iterative procedure for the solution of biharmonic equations. The method 
is also suitable for implementation on vector supercomputers. Notice that the APP and DPP 2 methods 
can also be applied with the generalized CG method (Axelsson 1980 and Eisentat et al. 1983) for 
non -symmetric matrix problems. Extension of the present method for general fourth order elliptic 
partial differential equations and practical problems of interest to computational fluid mechanics is 
under investigation; the results will be reported in a separate paper. 
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NUMBER Or TERMS 


Fig. 1 Comparison of two preconditionings for EQ = 9801 
(initial u° based on random numbers) 



C = 0.6282*1070.1905*10" = 3297.6 

Fig. 2 Eigenvalue spectrum of the original matrix B for EQ = 225 





0 1 0 20 30 40 0 10 20 30 40 

k=1, C = 0.3291*1070. 3265*10” = 1007.9 k=5. C = 0.3497*1070.2183*10° = 160.2 



0 10 20 30 40 0 10 20 30 40 

k= 10. C = 0.3565*1070.6638*10° = 53.7 k=15, C = 0.3669*1070.1368*10' = 26.8 


Fig. 3 Eigenvalue spectra of the preconditioned matrices BM' 1 for EQ = 225 
(DPP method for k = l, 5, 10, and 15) 
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k=10, C = 0.4935*1070.3884*10' = 12.7 k=15, C = 0.7147*1070.9526*10' = 7.5 


Fig. 4 Eigenvalue spectra of the preconditioned matrices BM' 1 for EQ=225 
(APP method for k = l, 5, 10, and 15) 
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